home *** CD-ROM | disk | FTP | other *** search
/ Pascal Super Library / Pascal Super Library (CW International)(1997).bin / MATH / NRPAS13 / FOUR1.PAS < prev    next >
Pascal/Delphi Source File  |  1991-04-29  |  2KB  |  57 lines
  1. PROCEDURE four1(VAR data: gldarray; nn,isign: integer);
  2. (* Programs using routine FOUR1 must define type
  3. TYPE
  4.    gldarray = ARRAY [1..nn2] OF real;
  5. in the calling routine, where nn2=nn+nn. *)
  6. VAR
  7.    ii,jj,n,mmax,m,j,istep,i: integer;
  8.    wtemp,wr,wpr,wpi,wi,theta: double;
  9.    tempr,tempi: real;
  10. BEGIN
  11.    n := 2*nn;
  12.    j := 1;
  13.    FOR ii := 1 TO nn DO BEGIN
  14.       i := 2*ii-1;
  15.       IF (j > i) THEN BEGIN
  16.          tempr := data[j];
  17.          tempi := data[j+1];
  18.          data[j] := data[i];
  19.          data[j+1] := data[i+1];
  20.          data[i] := tempr;
  21.          data[i+1] := tempi
  22.       END;
  23.       m := n DIV 2;
  24.       WHILE ((m >= 2) AND (j > m))  DO BEGIN
  25.          j := j-m;
  26.          m := m DIV 2
  27.       END;
  28.       j := j+m
  29.    END;
  30.    mmax := 2;
  31.    WHILE (n > mmax) DO BEGIN
  32.       istep := 2*mmax;
  33.       theta := 6.28318530717959/(isign*mmax);
  34.       wpr := -2.0*sqr(sin(0.5*theta));
  35.       wpi := sin(theta);
  36.       wr := 1.0;
  37.       wi := 0.0;
  38.       FOR ii := 1 TO (mmax DIV 2) DO BEGIN
  39.          m := 2*ii-1;
  40.          FOR jj := 0 TO ((n-m) DIV istep) DO BEGIN
  41.             i := m + jj*istep;
  42.             j := i+mmax;
  43.             tempr := sngl(wr)*data[j]-sngl(wi)*data[j+1];
  44.             tempi := sngl(wr)*data[j+1]+sngl(wi)*data[j];
  45.             data[j] := data[i]-tempr;
  46.             data[j+1] := data[i+1]-tempi;
  47.             data[i] := data[i]+tempr;
  48.             data[i+1] := data[i+1]+tempi
  49.          END;
  50.          wtemp := wr;
  51.          wr := wr*wpr-wi*wpi+wr;
  52.          wi := wi*wpr+wtemp*wpi+wi
  53.       END;
  54.       mmax := istep
  55.    END
  56. END;
  57.